The energetics of oxide surfaces by quantum Monte Carlo 
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Density functional theory is widely used in surface science, but gives poor accuracy for surface 
energetics in many cases. We propose a practical strategy for using quantum Monte Carlo techniques 
to correct DFT predictions, and we demonstrate the operation of this strategy for the formation 
energy of the MgO (001) surface and the adsorption energy of the H2O molecule on this surface. 
We note the possibility of applying the strategy to other surface problems that may be affected by 
large DFT errors. 



For many years, electronic-structure techniques have 
played a major role in surface and interface science. The 
most widely used of these techniques is density func- 
tional theory (DFT) [1], which has been employed to 
study many problems, including metal-oxide adhesion, 
catalysis and corrosion. Yet there is evidence that com- 
monly used DFT approximations are often seriously in 
error for basic quantities like surface formation energies 
and molecular adsorption energies, ft has been noted [2] 
that quantum Monte Carlo (QMC) techniques may be 
able to overcome these problems, because of their higher 
accuracy. We propose here a general strategy for using 
QMC to assess and correct DFT predictions for surface 
formation and molecular adsorption energies, based on 
the idea that the important DFT errors are localised near 
the surface, so that the comparison of QMC with DFT 
for small systems may often suffice to give the informa- 
tion needed. We will show the practical operation of this 
strategy for the formation energy of the MgO (001) sur- 
face and the adsorption energy of water on this surface, 
comparing where possible with experimental data. 

Insight into DFT errors for surface energetics comes 
from work on the jellium surface [3,4]. Jellium is the 
homogeneous interacting electron gas neutralised by a 
uniform background; its density n is characterised by the 
mean inter-electron distance r s , defined by (47rr 3 /3)n = 
1, with r s in atomic units. The planar jellium surface 
is formed by having the neutralising background occupy 
only the half-space x < 0, so that the electron number 
density n{x) in the ground state goes to its bulk value n 
for x — > — 00 and to for x — > 00. Accurate results for 
the formation energy of the jellium surface have been ob- 
tained [4] by extrapolating QMC calculations on neutral 
jellium spheres of different radius R. The extrapolation 
was performed by studying the large- R behaviour of the 
difference between the QMC total energy and the total 
energy calculated with DFT approximations, the main 
such approximations being: the local density approxima- 
tion (LDA) [1]; the generalised gradient approximation 
(GGA) in the Perdew-Burke-Ernzerhof form (PBE) [5]; 



and the meta-GGA [6]. This jellium surface work showed 
that: (i) for 2 < r s < 5, typical of simple metals such as 
Al and Na, the meta-GGA gives a very accurate surface 
energy a, followed closely by LDA, with GGA being too 
low by a significant amount; (ii) as r s falls below 2, the 
GGA errors rapidly worsen. In the region r s ~ 1.5, char- 
acteristic of transition metals and many oxides, the GGA 
a is too low by ~ 0.3 J m~ 2 , a serious error, because the 
surface energies of these materials are themselves in the 
region of 1 J m~ 2 . 

Indirect confirmation for large GGA errors in surface 
energies comes from a study [7] of the work of adhesion 
FU a dh of Pd (111) to CV-AI2O3 (0001), for which accurate 
measurements are available. (Here, Wadh is the reversible 
work per unit area needed to separate the system con- 
taining the oxide-metal interface into its metal and ox- 
ide consituents.) The GGA value W a dh — 1.6 J m~ 2 is 
far below the LDA and experimental values of 2.4 and 
2.8 J m~ 2 . The authors argue [7] that the GGA error 
comes mainly from errors in the free surface energies, 
and semiquantitatively relate these errors to GGA errors 
for the jellium surface. They suggest the general use of 
QMC jellium surface energies to correct DFT predictions 
for the surface energies of real materials [7] . 

The strategy we propose also uses QMC to correct 
DFT, but we apply QMC directly to the system of inter- 
est. In principle, QMC could be applied by brute force 
to the large slab systems commonly used to model sur- 
faces in DFT calculations. However, since QMC is far 
more costly than DFT, this is not generally feasible at 
present. It is also unnecessary, and not the best way of 
gaining insight. Since DFT errors for surface energet- 
ics are expected to be localised in the surface region [8], 
an accurate assessment of these errors should be given 
by QMC calculations only on the atoms in the surface 
region. This implies that thin slabs, containing only a 
few atomic layers, should suffice to assess the difference 
between the surface energy given by QMC and by DFT 
approximations, so that this difference will converge more 
rapidly with increasing slab thickness than the separate 
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surface energies. This strategy resembles that used to 
extract the energy of the jellium surface from calcula- 
tions on jellium spheres [4]. We propose to use the same 
scheme for molecular adsorption energies. The molecule 
is placed on the surface of a thin slab, and we study the 
difference between the QMC and DFT adsorption ener- 
gies, seeking convergence of this difference with increas- 
ing slab thickness. 

We have studied the practical feasibility of this strat- 
egy for the formation energy a of the MgO (001) 
surface. Our DFT calculations used the standard 
pseudopotential/plane- wave techniques [1], and were per- 
formed using the VASP code [9]. The surface was mod- 
elled using periodically repeated slab geometry, the cal- 
culation conditions being characterised by basis-set com- 
pleteness (plane- wave cut-off energy E cut ), Brillouin-zone 
sampling of the electronic states, the width L of the vac- 
uum layer separating successive slabs, and the number of 
layers Mayer in each slab. The surface formation energy 
is a = (E'siab - E hnlk )/A, with E sl&b the energy of the 
slab system, per repeating cell, -Ebuik the energy of the 
same number of atoms of the bulk material, and A the 
total surface area (both faces) of the slab, per repeating 
cell. This definition applies for all Mayor > 1- The bulk 
energy -Ebuik is Mayer times the bulk energy per layer 
ebuik, and it is convenient to obtain ebuik from the dif- 
ference of -Esiab values for successive values of Mayer in 
the limit of large Mayor- For given Mayer, we always in- 
sist on convergence of the calculated a with respect to 
_E cu t, BZ sampling and L to within 0.01 J m -2 (this tol- 
erance is satisfied for L > 6 A) . As expected from earlier 
work, ct for MgO (001) converges rapidly with respect 
to M ayer , the residual errors being below 0.01 J m~ 2 for 
Mayor > 2. For a given DFT approximation, the calcu- 
lated ct depends a little on MgO lattice parameter ao . For 
the experimental value ao = 4.21 A, we obtain a = 1.24 
and 0.87 J mT 2 with LDA and GGA(PBE) respectively. 
The difference of 0.37 J m~ 2 between the two is very sim- 
ilar to the difference of ~ 0.4 J m~ 2 between the LDA 
and GGA surface energies of a-Al 2 3 (0001) [7]. 

The calculation of ct by QMC is not standard, and 
we are not aware of previous calculations of ct for any 
oxide surface using QMC, though our recent QMC cal- 
culations on perfect and defective MgO crystals [10] in- 
dicated the feasibility of the present calculations. We 
refer the reader to reviews for details of QMC (e.g. [11]). 
We recall that for high-precision results it is essential 
to use diffusion Monte Carlo (DMC), in which the many- 
electron wavefunction is evolved in imaginary time, start- 
ing from an optimised trial wavefunction generated in 
prior variational Monte Carlo calculations. The only er- 
ror inherent in DMC is "fixed-node" error, due to the 
fact the nodes of the many-electron wavefunction arc 
constrained to be those of the trial wavefunction. For 
many systems, including jellium, the evidence is that 
fixed-node error is extremely small. For wide-gap sys- 



tems such as MgO, the errors should be no greater. Our 
calculations were performed with the CASINO code [12], 
using the same Hartree-Fock pseudopotentials as in our 
previous work [10]. The trial wavefunctions were of the 
usual Slater- Jastrow type, with single-electron orbitals 
obtained with the plane- wave code PWSCF [13], gener- 
ally using the large plane-wave cut-off of 4082 eV. These 
orbitals were represented in CASINO using the recently 
reported "blip-function" real-space basis set [14]. The 
DMC calculations all used a time-step of 0.005 a.u., and 
mean number of walkers equal to 10,240. The calcula- 
tions were done with free boundary conditions (i.e. no 
periodicity) normal to the surface. 

Our DMC calculations were performed on a series of 
MgO slabs with the number of layers Mayer running from 
1 to 5. For each Mayer, convergence must be demon- 
strated with respect to basis-set completeness and size of 
repeating surface unit cell. Basis-set errors with the blip 
basis set are readily made negligible, as shown earlier [14]. 
In DMC calculations, the wavefunctions are real, so that 
Brillouin-zone sampling is generally impossible, and cal- 
culations arc usually performed at the T-point; this is 
why convergence with respect to size of surface unit cell 
must be checked. Our main DMC calculations used the 
2x2 surface unit cell, for which the repeating cells con- 
tain from 16 (Mayer = 1) to 80 (Mayer = 5) ions; we show 
below that larger surface cells would give almost identical 
results. 

The raw output from these calculations is DMC to- 
tal energies 

E™c for the five iVlayor va i U es. Follow- 
ing our strategy, we now study the difference Ai? s i a b = 
E ?i^° - E EIb '. with thc DFT slab energy calculated with 
exactly the same slab and the same (r-point) BZ sam- 
pling as in the DMC calculations. Since the jellium re- 
sults indicate that LDA surface energies are likely to be 
closer to DMC than those from GGA, we use LDA values 
for E®^- . When plotted against Mayer, A-E s i a b will tend 
asymptotically (Mayer — * oo) to a straight line, whose 
slope is equal to the difference between the DMC and 
LDA bulk energies Aebuik per layer, and whose Mayer = 
intercept divided by A gives the difference of DMC and 
DFT surface energies Act ee ct dmc - ct dft . Since Ae bulk 
is large, and since Ai? s i a b contains the statistical errors 
of DMC, it is helpful to start this analysis by performing 
a least-squares straight-line fit a + bN\ ayel to the values 
of Ai? s iab, and then to use the resulting b value to form 
the quantity AE s ] ab = AE s i ab - b Mayor • The Mayor — > oo 
straight-line asymptote of AE s i a b has the same Mayer = 
intercept as that of Ai? s i a b- For Mayor = 1,2,... 5, 
we find the five values AE slah = -0.019(2), -0.009(6), 
-0.007(9), -0.011(13) and -0.014(15) J nT 2 . This im- 
mediately shows that the DMC and LDA values of ct 
arc almost exactly the same. Within our rather small 
statistical errors of at worst 0.015 J m~ 2 , thc difference 
between the DMC and LDA surface energies has the very 
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small value of —0.01 J m~ 2 . To check the errors due to 
use of the 2 x 2 surface cell (i.e. errors of BZ sampling), we 
have performed LDA calculations on slabs having large 
surface cells with a series of A^ aycr values, using T-point 
sampling. The a values thus obtained differ from the 
LDA a value extracted by similar T-point calculations 
with the 2 x 2 surface cell by only 0.01 J m -2 . The LDA 
value of a for the lattice parameter we are using, con- 
verged with respect to BZ sampling and slab thickness, 
is 1.20 J m -2 , and we conclude from our Ai? s i a b values 
that the fully converged DMC value is 1.19 ±0.01 J m~ 2 . 
MgO is one of the few oxides for which reasonably reli- 
able experimental values of the surface energy are avail- 
able [15]. Exploiting the fact that MgO cleaves readily 
along the (001) plane, the experiments measure the work 
of cleavage, thus ensuring that the results for a cannot 
be influenced by surface contamination. Our cxdmc value 
of 1.19 J m~ 2 is consistent with the measured values [15], 
which fall in the range 1.04 — 1.20 J m~ 2 . 

In applying QMC to correct DFT predictions for the 
adsorption energy E a d s of H 2 on MgO (001), we as- 
sume, in accord with experimental and theoretical indi- 
cations [16,17], that the molecule lies almost fiat on the 
surface, with the water O atom almost above a surface 
Mg ion, the water O-H bonds pointing towards surface 
O ions (Fig. 1). The adsorption energy is defined to be 
Eg.ds = £ , H 2 o + -Ebaro slab - -B s iab+H 2 o , where the terms on 
the right are the energy of the isolated H2O molecule, the 
energy of the bare MgO (001) slab, and the energy of the 
slab with the H 2 molecule adsorbed on the surface, all 
three systems being fully relaxed to equilibrium. In our 
DFT calculations, we require that E^ds be converged to 
within 10 meV with respect to plane- wave cut-off E cut , 
BZ sampling, and vacuum width L. Furthermore, since 
we want -E a ds for an isolated molecule, we examine the 
dependence of -E^ds on the size of the surface unit cell. 
We find that with the 2x2 cell, E a d s is already converged 
to better than 10 meV (this was tested by doing calcula- 
tions up to 5 x 5 surface unit cells) . With these tolerances 
always applied, we then study the dependence of i? a ds on 
the number of layers in the slab Mayer- As Mayer in- 
creases, i^ads ceases to change by more than 2 meV for 
Mayor > 2. Our calculated values of E a d s with LDA and 
GGA(PBE) are 0.92 and 0.43 eV respectively. 

In the DMC calculations, we obtained Eh 2 o using pe- 
riodically repeated cubes of different lengths d, with a 
single molecule in each cube, using the molecular geome- 
try taken from PBE (O-H bond length = 0.978 A, bond- 
angle = 104.4°). There is a weak dipole-dipole correction, 
going as d~ 3 , but we extrapolate to infinite d to obtain 
En 2 o with a technical error uncertainty of only a few 
meV. The slab energy -Ebare s i a b is taken from our DMC 
calculations on a (see above). For H 2 on the slab, strict 
application of our strategy would require us to use the 
relaxed configuration obtained from DMC calculations. 
We are not yet able to do this, since the calculation of 



ionic forces with DMC is problematic for the moment 
(though see Ref. [21]). Instead, we use a relaxed con- 
figuration from the DFT approximation which appears 
to reproduce DMC most closely. Our DFT calculations 
show that the relaxed height of the molecule above the 
surface differs by ~ 0.15 A between LDA and PBE. We 
calculated the DMC energy fo H 2 on the 1-layer slab, 
with a series of geometries on a linear path between the 
relaxed LDA and PBE geometries, and we find that the 
PBE geometry is very close to giving the lowest DMC 
energy. All our DMC results therefore refer to relaxed 
PBE geometries. 

Since DFT calculations of -E a ds are converged with re- 
spect to surface cell size for 2 x 2 cells, our DMC calcula- 
tions are all done with this surface cell. With A^ ayer = 1 , 
we find E^% c = 0.63(3) cV, E^ E = 0.48 cV, so that 
A£ ads ee E*% c - -Epg s E = 0.15 cV. For A layor = 2, the 
results are E^% c = 0.57(4) cV, £f;f E = 0.42 eV, so that 
A_E ads = 0.15 eV, identical to the Mayor = 1 value within 
statistical errors. Using the Mayer — > oo PBE value of 
0.43 eV, we thus estimate the Mayer — > oo DMC value as 
0.58(3) eV. 

For i? a ds, a comparison with experiment can only be 
indicative at present. Measurements of LEED isotherms 
and isobars for H 2 adsorption on MgO (001) as a func- 
tion of coverage, extrapolated to zero coverage gives 
E a ds = 0.52 ± 0.10 eV [18]. Temperature programmed 
desorption experiments show a peak at T = 235 — 260 K 
due to desorption of water at initial monolayer cover- 
age [16,19]. The standard Redhead analysis, using the 
commonly assumed frequency prefactor of 10 13 sec -1 , 
yields an effective _E a d s of 0.63 — 0.67 eV. However, this 
must include a significant contribution from attractive 
water-water interactions, so that E a ds for isolated H 2 
should be somewhat lower, and thus perhaps consistent 
with the LEED value of 0.52 ± 0.10 eV. To compare our 
DMC value with this, a correction for vibrational energies 
is needed. We have performed GGA(PBE) calculations 
on the adsorbed molecule, which indicate that adsorption 
lowers the symmetric and asymmetric stretch modes of 
H 2 by 12 and 9 THz respectively, and raises the bond- 
bending mode by 2 THz. The associated zero-point en- 
ergies raise the adsorption energy by 39 meV. Transla- 
tional and rotational energies of the H 2 molecule in 
free space, and vibrations and librations of the adsorbed 
molecule relative to the surface [20] raise and lower the 
adsorption energy by 64 meV and 167 meV respectively. 
Altogether, vibrational effects lower the adsorption en- 
ergy by 64 meV, so that our corrected DMC adsorption 
energy is 0.51 cV. Our calculated value thus appears to 
be consistent with the experimental evidence, but clearly 
a more elaborate analysis would be needed to make this 
comparison robust. Our comparisons make it fairly clear 
that LDA gives a serious overestimate of the adsorption 
energy, while GGA (PBE) is considerably more accurate. 

Our proposed strategy is therefore feasible and useful 
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for MgO surface energetics. Many other important prob- 
lems could be addressed in the same way, including the 
CK-AI2O3 (0001) surface energy mentioned earlier [7]; the 
computational effort needed for this case would similar to 
MgO. For this and other applications of the strategy, we 
believe that current progress in the calculation of ionic 
forces and structural relaxation with DMC [21] will be 
very helpful. 

The ability to validate QMC against experiment for 
molecular adsorption energies is limited by the lack of 
adequate techniques for putting modelling and data into 
close contact. Progress in the ab initio prediction of 
TPD spectra [22] is encouraging in this respect. Finally, 
we note the importance of quantum chemistry meth- 
ods. For wide-gap materials, techniques such as MP2 
and CCSD(T) should deliver an accuracy similar to that 
of DMC. Major improvements in the scaling of these tech- 
niques with number of atoms and size of basis sets [23] 
and to perform them in periodic boundary conditions [24] 
should make it possible to use them within our proposed 
strategy. 

In summary, we have described a practical general 
strategy for using QMC calculations to assess and cor- 
rect the errors of DFT approximations for the energetics 
of surfaces. Our calculations on the surface formation 
energy of MgO (001) and the adsorption energy of H 2 
on this surface confirm the feasibility and usefulness of 
the strategy. The results support earlier inferences from 
the energetics of the jellium surface that the GGA surface 
formation energy for this type of material is substantially 
too low, and that LDA is more accurate. However, for 
the molecular adsorption energy, the reverse is true, with 
LDA errors being much greater than those of GGA. 
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